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1.  INTRODUCTION 


Frequently  data  analysis  and  mathematical  modeling  require  compu¬ 
tation  of  derivatives  of  functions  given  at  discrete  points  only. 

For  instance,  temperature  dependent  heat  capacity  can  be  determined^  by 
differentiating  heat  content.  Dependence  of  the  rate  of  plant  photo- 

o 

synthesis  on  light  intensity  is  computed  by  measuring  the  amount  of 
absorbed  substances  and  by  differentiating  collected  data  numerically. 
Another  example  is  the  rate  of  elimination  of  a  drug  from  blood,  ob¬ 
tained  by  differentiating  the  measured  concentration  at  successive  time 
instances.  These  are  three  examples  from  rather  diverse  fields  of  in¬ 
quiry.  _  There  are  numerous  investigations  at  BRL  also  requiring  differ¬ 
entiation  of  tabular  data.  One  example  of  such  a  problem  is  computation 

of  target  speed  and  acceleration  from  measured  position  coordinates.* 2 * * 5 

Another  is  computing  a  drag  coefficient  from  the  differential  equations 
of  a  flight  trajectory  by  substituting  velocity  and  acceleration  values 
computed  by  differentiating  position  coordinates  derived  from  radar 
data.  Also  the  speed  and  acceleration  of  a  projectile  in  the  bore  of 

a. gun  have  been  determined5  by  computing  derivatives  of  discrete  posi¬ 
tion  data  obtained  by  a  microwave  interferometer.  Further  examples  of 
BRL  research  requiring  differentiation  of  tabular  data  include  deter¬ 
mination  of  the  slope  of  a  shock  front  from  discrete  position  points, 
computation  of  lead  angle  for  antitank  guns,  and  others. 

In  the  cases  mentioned  above  the  derivative  was  obtained  by  BRL 
researchers  either  by  drawing  a  smooth  curve  through  the  plotted  points 
and  then  measuring  its  slope,  or  by  a  moving  polynomial  arc  technique, 
or  sometimes  by  use  of  a  mechanical  spline. 

There  are  numerous  open  literature  publications  that  discuss 
various  numerical  methods  for  computing  derivatives  of  tabular  functions. 


T.F.  Dolgopolova  and  V.K.  Ivanov,  On  Numerical  Differentiation  Zh 
Vychisl  Mat.  Fiz.  6,  3,  570-576,  1966. 


2 

R.S.  Anderssen  and  P.  Bloomfield,  "A  Time  Series  Approach  to  Numerical 
Differentiation."  Technometrics  16,  1  69-75,  1974 


C.  Masaitis,  G.  Francis  and  V.  Woodward,  "Survival  vs.  Horsepower  per 
Ton  Test  Data  Analysis,"  BRL  MR  2518,  1975.  (AD  #000321 1 L)  H 


R.F.  Lieske  and  A.M.  Mackenzie,  "Determination  of  Aerodynamic  Drag  from 
Radar  Data,"  BRL  MR  2210,  1972.  y 


5 

F.  Yagi,  B.  Jansen,  L.B.  Kennedy,  W.C.  Taylor,  "Analysis  of  Interfero- 
Repon  1'n  th6  B°re  °f  Sn,a11  A™S-“  BRL 
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Some  of  these  simply  differentiate  polynomial  interpolation  formulas;^ 

others  use  least  squares  fits  of  the  data  by  trigonometric* * 7 *  or  alge- 

braic  polynomials.  Still  others  use  polynomial  splines.  This  last 
method  imposes  certain  smoothness  conditions  on  the  fitted  functions, 

and  similar  conditions  are  introduced  by  applying  Tikhonov's^  regu- 

,  .  ..  ,  1,11,12 

larization  procedure. 

Smoothness  conditions  of  the  spline  functions  and  of  the  regulari¬ 
zation  procedure  attempt  to  overcome  the  sensitivity  of  computed  de¬ 
rivatives  to  small  perturbation,  i.e.,  the  ill-posedness  of  the  dif¬ 
ferentiation  problem.  Because  of  this  sensitivity  simple  methods  such 
as  central  difference  or  moving  polynomial  arc  seldom  produce  satis¬ 
factory  values  of  derivatives  of  tabular  data  containing  measurement 
errors.  Under  these  circumstances  much  more  complicated  methods  are 
required,  such  as  those  of  regularization,  spline,  or  the  method  de¬ 
scribed  in  this  report. 

Since  tabular  data  contain  information  about  the  function  at  dis¬ 
crete  points  only,  and  a  derivative  is  the  limit  of  a  certain  ratio,  it 
cannot  be  obtained  from  numerical  data  without  additional  assumptions. 
The  most  common  assumption  is  that  tabular  data  are  values  of  a  certain 
function  which  can  be  identified  on  the  basis  of  the  data  and  subse¬ 
quently  differentiated  to  yield  the  values  of  the  derivative.  For  in¬ 
stance,  the  approximation  of  a  derivative  by  the  divided  difference 
assumes  that  tabular  data  can  be  represented  by  a  straight  line,  at 
least  locally.  An  extension  of  this  assumption  is  obtained  by  fitting 


^D.B.  Hunter,  "An  Iterative  Method  of  Numerical  Differentiation,"  Comp. 

J.  3,  270-271,  1960. 

7A.  Talmi  and  G.  Gilat,  "Method  for  Smooth  Approximation  of  Data," 
Journal  of  Comp.  Physics  23,  93-123,  1977. 

^H.C.  Hershey,  J.L.  Zakin  and  R.  Simha,  "Numerical  Differentiation  of 
Equally  Spaced  and  Not  Equally  Spaced  Experimental  Data,"  Ind.  Eng. 
Chem.  Fundam.  6,  413-421,  1967. 

9P.  Craven  and  G.  Wahba,  "Smoothing  Noisy  Data  with  Spline  Functions," 
Numer.  Math.  31,  377-403,  1979. 

^9A.N.  Tikhonov,  "Solution  of  Incorrectly  Formulated  Problems  and  the 
Regularization  Method,"  Soviet  Math.  Dokl .  4,  1035-38,  1963. 

^J.  Cullum,  "Numerical  Differentiation  and  Regularization,"  Siam  J. 
Numer.  Anal .  8,  254-265,  1971 . 

1  2 

R.S.  Anderssen  and  P.  Bloomfield,  "Numerical  Differentiation  Pro¬ 
cedures  for  Non-Exact  Data,"  Numer.  Math.  22,  157-182,  1974. 
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a  polynomial  of  prescribed  degree  and  then  differentiating  it.  For 

computational  convenience,  orthogonal  polynomials  may  be  selected8 
and  data  represented  by  a  linear  combination  of  these  polynomials.  The 
moving  polynomial  arc  procedure  assumes  that  the  data  can  be  locally 
represented _ by  a  polynomial  of  a  specified  degree.  Thus,  this  pro-' 
cedure  is,  in  a  certain  sense,  a  generalization  of  the  divided  difference 
approximation,  which  represents  the  data  by  a  first  degree  polynomial 
fitted  to  two  data  points.  Representation  of  data  by  spline  functions, 
which  are  polynomials  on  the  intervals  between  the  nodes,  also  belongs 
to  the  same  class,  with  a  special  data  fitting  criterion  assuring  that 
the  L,,  norm  of  a  derivative  of  certain  order  is  sufficiently  small. 

Besides  polynomials  and  splines,  frequently  trigonometric  polyno¬ 
mials  and  linear  combinations  of  exponentials  (exponential  polynomials) 
are  used  to  approximate  tabular  data,  and  then  their  derivatives  are 
computed  by  differentiating  the  approximating  polynomials.  Also, 
approximation  of  tabular  data  by  a  linear  combination  of  arbitrary 

orthogonal  functions  is  frequently  discussed  in  the  literature.7  How- 
ru6[’  !n  aPP1ications  these  functions  are  nothing  more  than  Legendre  or 
Chebychev  polynomials,  trigonometric  functions,  or  just  powers  of  the 
independent  variable  that  are  not  even  mutually  orthogonal. 


iriub  almost  an  tne  approximations  encountered  in  practice  are 
least  squares  fits  by  algebraic,  trigonometric,  or  exponential  polyno- 
mla  ..u  a  Qlobal  type  or  local  as  with  a  moving  polynomial  arc, 

and  either  without  additional  constraints  or  with  certain  constraints 
as  in  the  case  of  the  spline  procedure.  Consequently,  all  these  methods 
can  be  generalized  by  choosing  an  approximating  function  from  the  al¬ 
gebra,  call  it  A  for  convenience,  generated  by  algebraic,  trigonometric, 
and  exponential  polynomials  and  by  properly  formulating  an  approximation 
criterion. 


Usually  the  type  of  linear  combination  of  selected  functions  is 
chosen  in  advance  and  only  the  weighting  factors  of  this  combination 
are  adjusted  to  the  data  according  to  a  selected  criterion  (mostly  the 
least  sum  of  residuals  squared).  Such  a  choice  presupposes  that  the 
type  of  the  function  represented  by  the  tabular  data  is  known  independ¬ 
ently  of  the  data,  and  the  latter  are  used  only  to  determine  the  weight¬ 
ing  factors.  A  rather  rare  deviation  from  this  procedure  compares  the 
reduction  of  the  root  mean  squared  error  (RMSE)  of  the  fits  by  polynomials 


of  increasing  degree  and  selects  an  approximating  polynomial  accordingly, 
In  this  report  a  procedure  is  described  for  automatically  selecting  an 
element  of  the  algebra  A  that  minimizes  the  RMSE  in  an  autoregressive 
model  subject  to  constraints  implied  by  three  additional  assumptions. 


All  the  approximations  mentioned  above  yield  functional  values 
y(t)  as  linear  combinations  of  tabular  data,  x.'s,  with  the  weighting 

factors  being  linear  combinations  of  the  selected  basis  functions  such 
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as  polynomials,  exponentials,  etc.,  say  4>-|(t),  cj>2  ( t ) , . . .  t) : 

y(t)  =  I  x  l  C  <f>.(t).  (1.1) 

j  J  i=l  J1  1 


The  coefficients  C..  may  be  independent  of  the  data  as  in  moving  poly- 

J 

nomial  arc  or  regression  procedures,  or  they  may  be  rather  complex 
functions  of  data  points  as  in  regularization  and  spline  methods.  Dif¬ 
ferentiation  of  (1.1)  yields  the  derivative  of  the  approximating  func¬ 
tion  as  a  linear  combination  of  data  values: 

/  (t)  =  I  x.  I  C..  <J>',(t).  (1.2) 

j  J  i  J 

The  weighting  factors  in  (1.2)  depend  on  t,  i.e.,  on  the  position  of 
the  argument  t  relative  to  the  data  points,  t/s,  used  to  determine  the 

coefficients  C...  If  the  position  of  t  is  fixed,  say,  at  the  midpoint 

J  * 

of  the  span  of  t/s,  then  the  coefficients  C..  are  independent  of  t. 
j  j  i 

However,  in  either  case  the  derivative  is  expressed  as  a  linear  combina¬ 
tion  of  the  data  values,  just  as  in  the  special  case  of  its  approximation 
by  a  divided  difference. 

The  regularization  procedure  ’  *  obtains  a  numerical  derivative 

by  selecting  an  element  in  a  Sobolev  space  by  minimizing  the  sum  of 
residuals  squared  plus  a  term  proportional  to  the  square  of  the  norm  of 
the  approximating  function.  Thus,  this  is  an  extension  of  the  spline 

q 

approximation  which  considers  only  a  particular  norm  of  Sobolev  space. 
Furthermore,  in  practical  applications  the  regularization  procedure  is 
restricted  to  readily  computable  elements  of  the  Sobolev  space  and  typi- 

11  12 

cally  considers  only  trigonometric  polynomials.  ’ 

In  view  of  this  it  is  natural  to  consider  only  representations  of 
derivatives  of  tabular  functions  by  linear  combinations  of  the  data 
values  and  to  determine  the  weighting  factors  of  such  a  combination 
directly  from  the  data  with  a  minimal  set  of  assumptions.  This  is  the 
approach  of  the  present  report. 

Our  basic  assumption  is  that  the  given  data  represent  an  element 
of  the  algebra  A  as  defined  above.  A  criterion  as  to  what  constitutes  the 
best  representation  of  the  data  is  our  second  assumption.  It  depends  on 
the  sum  of  the  squares  of  certain  residuals  and  on  the  degree  of  ill-con¬ 
ditioning  of  the  problem.  The  definition  of  this  criterion  is  contained 
in  the  next  section.  The  criterion  yields  only  the  structure  of  an  ap¬ 
proximating  element.  The  values  of  the  element  parameters  can  be  subse¬ 
quently  determined  by  a  linear  least  squares  fit  to  a  selected  segment 
of  the  data.  A  derivative  of  any  order  is  obtained  either  by  differ- 
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entiating  such  a  fitted  approximation  or  by  computing  the  weighting 
factors  directly  from  the  structure  of  the  element. 

Our  third  assumption  is  on  the  number  of  data  values  (and  weight¬ 
ing  factors)  in  the  linear  representation  of  a  derivative.  It  is  shown 
in  this  report  that  the  derivative  of  every  element  of  A  is  representa¬ 
ble  as  a  linear  combination  of  its  values  and,  conversely,  a  function 
is  in  A  if  its  derivative  can  be  represented  as  a  linear  combination  of 
a  finite  number  of  its  equally  spaced  values  from  any  of  the  several 
data  intervals. 

Our  fourth  assumption  is  that  the  step  size  in  the  tabular  data  is 
sufficiently  small  so  that  the  structure  of  the  corresponding  element 
of  A  is  uniquely  defined  by  the  values  of  the  data.  This  assumption 
implies  that  the  selected  element  of  A  is  sufficiently  smooth,  with  the 
degree  of  smoothness  determined  by  the  step  size  and  the  accuracy  of 
the  data. 

These  four  assumptions,  spelled  out  precisely  in  the  following 
sections,  yield  values  of  the  derivatives  that  are  sufficiently  faithful 
to  the  empirical  data  and  at  the  same  time  provide  good  approximations 
of  derivatives  when  applied  to  synthetic  data.  Sections  2  and  3,  to¬ 
gether  with  Appendix  A,  present  formal  derivations  of  approximating  for- 
ulas  for  both  tabular  data  and  related  derivatives,  first  determining 
suitable  structure  based  on  eigenvalues  and  then  finding  parameters  by 
least  squares  procedures.  Section  4  deals  with  the  frequently  encount¬ 
ered  case  in  which  no  eigenvalue  is  repeated.  Here  the  formulas  for  the 
weighting  coefficients  can  be  simplified  and  rewritten  in  closed  form  in 
terms  of  real  elements  only,  with  a  consequent  reduction  in  round-off 
error  as  well  as  enhanced  practical  convenience.  Section  5  presents 
formulas  for  approximating  theoretical  error  bounds  for  derivatives 
computed  by  the  procedures  of  Section  4.  Appendix  B  illustrates  the 
accuracy  of  the  method  and  compares  it  to  the  accuracy  of  derivatives 
obtained  by  other  commonly  used  methods. 


2.  APPROXIMATION  OF  TABULAR  FUNCTIONS 
AND  THEIR  DERIVATIVES 

As  stated  in  the  introduction  we  approximate  tabular  data  by  an 
element  of  the  algebra  A  generated  by  algebraic,  trigonometric,  and 
exponential  polynomials  of  a  real  variable  t  e[0,T]  with  real  coeffi¬ 
cients.  An  approximating  element  of  A  is  selected  in  two  steps,  by  first 
determining  its  structure,  most  appropriate  in  a  certain  sense,  and  next 
by  determining  its  parameters  by  a  least  squares  fit. 

2.1  Data  Structure 


The  structure  of  an  element  of  A  approximating  tabular  data  is 
obtained  by  observing  that  a  function  which  satisfies  a  family  of  linear 
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difference  equations  dependent  on  the  scaling  of  the  independent  varia¬ 
ble  is  contained  in  A  and,  conversely,  every  element  of  A  satisfies  a 
linear  difference  equation  of  a  fixed  order  with  the  coefficients  de¬ 
pendent  on  the  scaling  of  the  independent  variable. 


By  Demoivre's  theorem  trigonometric  polynomials  can  be  expressed 

t  t  • 

as  linear  combinations  of  expressions  (cos  9.  ±  i  sin  0.)  =  Xj  with 

±i0.  ... 

X.  =  e  J  ,  i.e.,  a  trigonometric  polynomial  is  a  linear  combination 

of  exponentials  with  complex  bases.  Similarly,  exponential  polynomials 


e,t 


are  of  the  form  Y  c .  e  J  =  Y  c.  X..  Consequently,  every  element  of  A 
L  J  JO 

is  representable  in  the  form 


m 


"j 


f(t)  -  l  l 

j=l  p=0 


c.  Xt  tP 
OP  J 


(2.1) 


with  X/s  either  complex  or  real.  To  every  complex  X. 

J  .  .  J 

corresponds  its  complex  conjugate,  say,  X^  with  n^ 


in  (2.1)  there 

=  n ..  Let 
J 


A  >  0  and  write  (2.1)  in  the  form 

m  nj  t  P 

d  (- 

_  ^ iD  VA' 

0=1  P; 


m  J  f  K  a  A 

fit)  -  I  I  c'  (|)  , 

j=l  p=0  JP  A  J 


(2.2) 


where  d.  =  c.  Ap.  Define  the  polynomial  PA(X)  of  degree  k  where 


k=  Y  (n.+l)  (2-3) 

0=1  0 

as  follows 

m  .  n .+1 

P  (X)  =  n  (X-X,A)  J  .  (2.4) 

A  0=1  J 

Let  the  operator  B,  be  defined  by  B^  f(t)  =  f(t-A).  Then  it  follows 

from  the  properties  of  linear  difference  equations  with  constant  co¬ 
efficients  that  f(t)  given  by  (2.2)  satisfies  the  k-th  order  difference 
equation 

PA(BA)f(t)  =  0,  (2.5) 

i.e.,  if  f(t)eA  then  for  every  A  >  0  f(t)  satisfies  the  difference 
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equation  of  the  form  (2.5)  defined  by  the  A^,  their  multiplicities 
nj+l  >  and  A.  This  can  also  be  shown  by  substituting  (2.1)  in  (2.5). 

The  converse  is  also  true,  i.e.  if 


(a)  y(t)  satisfies  (2.5)  for  every  A  >  0  and  some  A/s  and 

J 

(b)  y(t)  has  bounded  derivatives  up  to  order  k  on  the 
interval  [0,T] 

then  y(t)eA.  Without  loss  of  generality  we  can  assume  that 


(°)  y(t)  satisfies  no  linear  difference  equation  of  order 
less  than  the  k  given  by  (2.3).  The  fact  that  such 
a  y(t)  is  in  A  follows  from  the  following  proposition: 


Proposition  1.  If  y(t)  satisfies  conditions  (a),  (b),  and 

(c),  then  it  is  a  solution  of  the  differential  equation 


P(D)y(t)  =  0 


(2.6) 


where  D  =  and 


m  n  .+1 

P(K)  =  n  U-  logA .)  J 
J-l  J 


(2.7) 


0-b A)  y 


Proof:  By  writing  y^  =  lim  '  A 


A-*0  A 


j; -  and  by  applying  Lemma  1  in 


Appendix  A  we  have  y^  =  jl  ,^N(A) 

k‘  dAk 


A=0 


.  Hence  by  (A. 10)  and 


Lemma  3  of  Appendix  A 


(t>  -  i  (-i)J-1 
j-i 


(t) . 


Consequently,  by  the  definition  of  S.  .  (see  Lemma  3),  y(t)  satisfies 
/  _  _  x  ^  »  J 


In  view  of  the  relation  between  the  algebra  A  and  the  equation  (2.5) 
an  element  of  A  that  approximates  tabular  data  can  be  selected  by  con¬ 
structing  an  appropriate  autoregressive  model  by,  say,  a  procedure  de¬ 
veloped  for  time  series  analysis^.  However,  this  requires  an  assump- 


G.E.  Box  and  G.M.  Jenkins,  Time  Series  Analysis:  Forecasting  and 
Control.  Revised  Edition.  Holden-Day,  1976. 
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tion  of  stationarity  or  else  instead  of  the  original  data,  say  x(n), 

n=0,l,...,N,  one  must  approximate  differences  A^x(n)  of  order  d  that  are 
assumed  to  be  stationary.  However,  this  assumption  need  not  hold  for 
tabular  data  that  must  be  differentiated.  For  instance,  if  the  under¬ 
lying  function  is  exponential  the  differences  of  any  order  are  exponential 
and  hence  non-stationary.  An  example  of  this  type  of  data  is  a  concen¬ 
tration  of  injected  drug  in  blood  measured  as  a  function  of  time  after 
injection.  In  view  of  this,  instead  of  attempting  to  determine  the 
order  of  differences  that  may  produce  stationary  series  and  at  the  same 
time  considering  possible  periodicity  ("seasonal"  variation),  we  choose 
a  direct  method  for  estimating  the  coefficients  of  an  autoregressive 
model.  To  this  end  we  write  x(n,p,q)  for  x(p+qn)  where  q  is  any  positive 
integer  and  p=0,l,2,... ,q-l,n=0,l,... ,Np,  where  Np=[(N-p)/qJ  is  the  largest 

integer  not  exceeding  (N-p)/q.  Let  y(t)eA  be  an  approximating  function 
of  the  data ,  i .e. 

x(n ,p,q)  =  y(rA)  +  er,  (2.8) 

where  r=p+qn  and  e  is  an  observation  error  assumed  to  be  weakly  sta- 

r  2 

tionary  white  noise  with  zero  mean  and  variance  a  .  The  function  y(t) 
satisfies  (2.5)  for  a  suitable  polynomial  P,  say,  of  degree  k  since 
y(t)eA.  We  write  this  equation  as  follows 

k 

y(rA)  =  l  a.  y[(r-jq)A]  ,  (2.9) 

j=l  3 


where  the  a.'s  remain  to  be  determined.  By  substituting 
we  get 

k 


x(n,p,q)  -  e  =  l  a,  [x(n-j,p,q)  -  e  .] . 

r  j=1  j  r-j 


(2.8)  in  (2.9) 
(2.10) 


By  transposing  the  terms  in  (2.10)  and  squaring  both  sides  we  get: 


1  x(n-j,p,q)  2 
j=l  J  J 


=  e 


k 

I 

j=l 


2  2 
3j  er-j 


(2.11) 


where  is  a  linear  combination  of  products  ee  with  u  f  v.  Since 
by  assumption  is  white  noise,  we  have  E(Pr)=0.  Thus,  by  taking 
expected  values  of  both  sides  of  (2.11)  we  get: 
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k 

2 

E 

x(n,p,q)  -  l  a.  x(n-j,p,q) 

L  J=1  J  J 

=  a2  +  a2  l  a2  (2.12) 

j=l  J 


We  replace  the  expected  value  of  the  left  hand  side  by  its  estimate 
(average)  and  get: 


N  p=0 


J- 

n=k+l 


x(n,p,q) 


l  a,  x(n-j ,p,q) 
j=l  J 


2  2  2 

=  G  +  <*  l  a,  ,  (2.13) 

j=l  J 

~  q;] 

where  N  -  £  (N  -  k)  .  Thus  we  get  from  (2.13): 

p=0  p 


o  = 


1  q-1  N 

t  lP 

IN  p=0  n=k+1 


x(n,p,q)  -  l 
: j~l 


x(n-j,p,q)l 


k 

1  +  -Ji 
J=1 


(2.14) 


We  choose  the  estimates  of  the  a.'s  that  minimize  the  variance  a2. 

J 

These  estimates  are  obtained  by  the  following  iterative  procedure.  Let 
a  denote  the  vector  (a^).  Assume  that  the  value  of  a,  say  a^^has 
been  obtained  on  the  u-th  iteration.  By  substituting  this  value  in 

r\ 

(2.14)  we  obtain  an  estimate  of  of  the  variance  of  e  .  Substitution  of 

ip  u 

this  value  for  0  in  the  second  term  of  the  right  hand  side  of  (2.13) 
yields: 


a 


2 


l  l  x(n,p,q)  -  l  a.  x(n-j ,p,q) 
p=0  n=k+l  L  j=i  J 


(2.15) 
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The  next  estimate  a^u+  ^  of  the  vector  a  is  obtained  by  minimizing  the 
right  hand  side  of  (2.15).  Partial  derivatives  of  this^expression  with 

respect  to  the  a.'s  when  set  to  zero  and  multiplied  by  N/2  yield 
J 

,  i  N 

k  q-1  p  ~  2 

l  a.  I  l  x(n-j,p,q)x(n-i,p,q)  -  N  a  a. 
j=l  J  p=0  n=k+l 


p=0 


np 

l  x(n,p,q)x(n-i ,p,q) 
n=k+l 


(2.16) 


i=l,2,...,k.  Thus,  if  M  is  the  matrix  of  the  normal  equations  of  the 
overdetermined  system 


l  a.  x(n- j ,p,q)  =  x(n,p,q),  (2.17) 

j=l  J 


n=k+l,...,N  ,  p=0,l ,. . .  ,q-l ,  and  X  is  the  right  hand  side  vector  of 

P  (u+1)  ...  2  . 

these  normal  equations,  then  the  vector  av  that  minimizes  a  in 

(2.15)  satisfies  the  equation 


(M  -  N  oj-  I)a(u+1  ^  =  X, 


(2.18) 


where  I  is  the  k  X  k  identity  matrix.  Successive  iterates  of  the  vector  a 

are  obtained  by  substituting  oZ  =  0  in  (2.18)  and  then  by  iterating 

(2.14)  and  (2.18).  Thus,  we  need  to  obtain  the  normal  equations  of  the 
overdetermined  system  (2.17)  and  then  solve  the  linear  equations  (2.18). 

In  order  to  compute  the  a.'s  by  this  procedure  we  have  to  choose 

J 

k  and  the  integer  q  in  (2.17).  Obviously  a  larger  number  of  model 
parameters,  i.e.  larger  k,  yields  a  model  better  matching  the  data,  A 
smaller  value  of  q  provides  a  larger  number  of  data  points,  i.e. 
describes  the  structure  of  the  data  more  accurately.  However,  in¬ 
creasing  k  as  well  as  reducing  q  makes  the  system  (2.18)  ill  conditioned. 
Hence  the  value  of  k  and  the  data  spacing  parameter  q  must  be  chosen  to 

minimize  oZ  in  (2.14)  and  at  the  same  time  to  prevent  the  matrix  in 
(2.18)  from  becoming  nearly  singular.  Thus,  we  have  two  conflicting 
criteria  for  selecting*  the  optimal  pair  (k,q).  As  usual,  a  measure 

n  _  . 

This  is  similar  to  the  solution  of  the  numerical  differentiation  prob¬ 
lem  by  regularization^  where  increasing  the  regularization  parameter 
reduces  ill-conditioning  of  the  problem  and  decreasing  the  parameter 
yields  a  better  fit  of  the  data. 
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of  optimality  must  be  chosen  heuristical ly.  Our  choice  is  an  index 

J(M)  =  a2(k,q)/[D(k,q)]^  ,  (2.19) 

where  D(k,q)  is  the  absolute  value  of  the  determinant  of  the  last  itera¬ 
tion  of  (2.18)  corresponding  to  the  choice  of  k  and  q.  Thus,  we  compute 
the  a^.'s  and  J(k,q)  for  k=l,2,...,kg  and  q=l,2,...,qg  and  select  the 

pair  (k,q)  and  the  corresponding  a/s  that  minimize  J(k,q).  We  impose 

an  additional  constraint  on  (k,q)  in  order  to  prevent  a  choice  of  a 
model  for  which  the  data  are  inadequate,  i.e.  a  model  that  contains 
terms  of  higher  frequency  than  can  be  determined  by  the  frequency  of 
the  data  points.  Thus,  if  wm  is  the  maximum  frequency  of  the  selected 

model  then  we  must  at  least  have 


q 


A  < 


2tt 

w 

m 


(2.20) 


Suppose  further  that  for  some  q  the  coefficients,  a/s,  in  (2.17)  yield 

a  real  negative  eigenvalue,  say  A.  <  0.  Then  the  term  c ..  Xn.  in  (2.1) 

.  i  »n  J  JO  J 

is  equal  to  c^  | j  cos  nir  for  every  n.  The  frequency  of  this  term 

is  tt  radians  per  qA  sec  or  ir/qA  radians  per  sec,  i.e.  we  have  w  >rr/qA, 

m 

contrary  to  requirement  (2.20). 


If  for  some  k  and  q  the  equation  (2.17)  has  an  eigenvalue  with  a 
negative  real  part,  say  A.  =  -a+ib  (a>0),  then  the  corresponding  term 

cj0  i°  (2.1)  is  expressible  as  c^.g  exp(pn+iwn)  where  cos  w  =  -a/Va2+b2, 

i.e.  w  >  j  if  expressed  in  radians  per  unit  time  equal  to  qA  sec.  There¬ 
fore  this  choice  of  q  yields  a  spacing  qA  with  less  than  four  data  points 
per  period  of  the  corresponding  term  in  (2.1).  Although  theoretically 
two  points  per  period  may  be  adequate  to  determine  the  real  and  imag¬ 
inary  parts  of  the  corresponding  eigenvalue,  even  three  points  per 
period  are  inadequate  when  the  data  contain  measuring  errors.  Further¬ 
more,  a  negative  real  part  only  implies  that  the  corresponding  frequency 
is  greater  than  tt/2  per  unit  time.  It  may  also  be  greater  than  it  and 
less  than  3tt/2,  in  which  case  the  spacing  qA  provides  less  than  two 
points  per  period.  This  is  the  reason  why  the  pairs  (k,q)  leading  to 
complex  roots  with  negative  real  part  are  rejected. 

In  summary,  the  models  (2.17)  are  determined  for  q=l ,2 . q„, 

k=l,2,...,kg  and  among  those  with  eigenvalues  having  non-negative  real 
components  that  one  which  yields  minimum  J(k,q)  in  (2.18)  is  selected. 
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When  the  data  is  very  noisy  this  selection  of  q  may  lead  to  a 
rather  large  step  size  qA  and,  thus,  may  eliminate  high  frequency  terms 
present  in  the  data  even  if  the  original  spacing  A  is  adequate  to  repre¬ 
sent  this  high  frequency.  This  may  happen  when  the  amplitudes  of  high 

frequency  terms  are  too  small  relative  to  the  measuring  error  e  to 

r 

be  determined  by  the  data  taken  at  any  spacing  A.  The  procedure  de¬ 
scribed  above  is  intended  to  determine  only  the  terms  of  (2.1)  for 
which  both  the  spacing  and  also  the  accuracy  of  the  data  are  adequate, 
and  this  works  satisfactorily  in  practice. 

2 . 2  Smoothing  and  Differentiation 


The  procedure  just  described  determines  the  coefficients  aj  in 

equation  (2.9),  i.e.  determines  the  structure  of  an  element  of  A  that 
minimizes  J(k,q).  The  corresponding  element  of  A  is  then  given  by  (2.1), 


i.e. 

of 


by  y(tQ  +  nqA) 


m 


I 

j=l 


y 


k 


l  c..  n1  y1?,  where  the  y.'s  are  the  roots 
i  =0  J1  J  J 


(2.21) 


of  multiplicity  n.+l.  By  reverting  back  from  qA  to  the  original  spacing 

J 

A  of  the  data  and  by  writing  X.  for  yV^  we  can  write  (2.21)  as  follows: 

J  J 


nj 


I  l 

j=l  i=0 


xu  » 


(2.22) 


where  c..  =  c^./q1 ,  u  =  |-  ,  and  y(tQ+t)  is  replaced  by  the  data  value 
xu=x(tQ/A  +  t/A). 


The  coefficients  in  (2.22)  are  selected  to  minimize  the  RMSE  of  the 
resulting  approximation.  This  may  be  either  a  global  or  a  local  approxi¬ 
mation  over  the  span,  say,  from  u  =  -K  to  u  =  K.  In  the  latter  case 
the  coefficients  c..  are  obtained  by  the  least  squares  method  from 

\J  * 

(2.22)  with  u  =  -K,  -K+1,..,K.  The  corresponding  normal  equations  are 


K  m  nj 

i  i  i 

u=-K  j=l  1=0 


„  ,,i+a  ,u  .u 

cji  u  xi  Ae  * 


U=-K 


*u  ““ 


or  in  matrix  notation 
S  C  =  Z. 


(2.23) 
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m 


C  =  S-1Z.  (2.24) 

Let  A  be  a  column  vector  with  components  v1^.  Then  substitution  of  C 
from  (2.24)  in  (2.22)  yields  a  smoothed  value  of  x 


(2.24) 


v 


Let  be  the  column  vector  with  components  uaAg.  Then  replacing  Z  in 

(2.25)  we  get: 


(2.26) 


u  and  v  but  not  on  x  ,  i.e. 


(2.26)  provides  a  smoothed  value  of  the  data  point  xy  as  a  linear  com¬ 
bination  of  the  data  values  xy.  If  (2.26)  is  used  to  compute  the 
smoothed  value  at  the  midpoint  of  the  span  of  the  data  from  -K  to  +K, 
(i.e.  xQ,  then  the  weighting  factors  remain  constant  as  the  span,  to¬ 
gether  with  its  midpoint,  is  translated  along  the  axis  of  the  independ¬ 
ent  variable  t. 


The  vector  C  given  by  (2.24)  can  be  used  to  compute  derivatives  of 
tabular  data.  By  differentiating  (2.1),  rescaling  the  independent 
variable  t  as  in  (2.22),  substituting  C  from  (2.24),  and  interchanging 
matrix  multiplication  and  summation  we  obtain  the  following  expression 
for  the  derivative: 


(2.27) 
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Similarly,  higher  derivatives  are  obtained  by  the  formula 


where  A 


(2.28) 


Weighting  factors  in  a  formula  like  (2.28)  for  derivatives  of 
various  orders  can  be  obtained  directly  from  the  structure  of  the  data 
defined  by  the  eigenvalues  of  a  corresponding  autoregressive  model 
without  computing  the  vector  C  in  (2.24).  Such  a  procedure  is  espe¬ 
cially  convenient  when  all  the  eigenvalues  are  simple,  and  this  case  is 
described  in  Section  4. 


3.  DERIVATIVES  AS  LINEAR  COMBINATIONS  OF 
FUNCTIONAL  VALUES 


It  was  mentioned  in  the  introduction  that  all  the  common  methods 
for  numerical  differentiation  express  derivatives  as  linear  combinations 
of  discrete  functional  values  in  the  form 


k 

/  (t)  =  l  a,y[t+(p-i )A] 
i  =  l  1 


(3.1) 


where  p  is  some  integer  and  A  is  the  spacing  between  data  points  used. 
Various  methods  apply  different  criteria  by  which  the  coefficients  a. 

in  (3.1)  are  selected.  A  generalization  of  all  these  methods  could  be 
obtained  by  identifying  a  set,  say  A,  of  functions  that  satisfy  the 
relation  such  as  (3.1)  and  by  defining  a  set  of  criteria  for  selecting 
the  coefficients.  This  we  do  below  and  in  Section  4. 

Obviously,  any  function  that  has  a  derivative  at  at  least  one  point 
satisfies  (3.1)  for  at  least  one  value  of_t.  Therefore  we  restrict  the 
set  A  by  insisting  that  every  element  of  A  satisfies  (3.1)  for  all  t's 
in  an  interval  such  as  [0,T].  If  we  allow  the  coefficients  a.,  to  vary 

with  t  in  a  different  fashion  for  every  function  in  A,_then  every  dif¬ 
ferentiable  and  non-vanishing  function  in  [0,T]  is  in  A,  since  we  can 
write  f'  (t)  =  [  f  (t)/f (t)] f (£)  which  is  a  special  case  of  (3.1). 
Therefore  we  define  the  set  A  as  a  set  of  functions  y(t)  each  of  which 
satisfies  (3.1)  with  a.'s  dependent  only  on  y(t),  p,  and  A  but  not  on 
t.  1 

We  choose  in  (3.1)  p=k+s+l  and  replace  the  summation  index  i  by 
-j+k+s+1 .  This  yields 
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(3.2) 


k+s 

y'  (t)  =  I  a(j,s)  y(t  +  jA), 
j=s+l 

where  a(j,s)  =  a_j+k+s+1 .  Thus,  the  set  A  is  defined  as  a  set  of 

IMi  -tl  ia*isiVkA2)  ,or  every  4  e[0'T]  and  every  1nteger  s 

The  following  proposition  is  implied  by  Lemmas  4  and  5  in  Appendix 


Proposition  2. 
braic,  trigonometric, 
[0,T]. 


A  =  A,  where  A  is  the  algebra  generated  by  alge- 
and  exponential  polynomials  over  the  interval 


and°write1rSt  ^  Sh°W  that  ^  C  A'  To  this  end  we  choose  s»-l  in  (3.2) 

yx(t)=  l  b.^)  y[t+(j-l)A]  ,  (3.3) 

j=l  J 

where  b^.'  =  a  ( j  - 1 ,  -1).  We  observe  that 

y(q,(t)  .  j  (3.4) 

J 


for  q  <  k  and  some  constants  b^q).  This  relation  is 

q-  Indeed,  (3.4)  holds  for  q=l  because  of  (3 
(3.4)  holds  for  q  <  k.  Differentiation  of  (3.4)  with 
yields: 


proved  by  induc- 

3).  Suppose 
respect  to  t 


(t) 


I  b,^  y'[  t+(j-l ) A] 

j=l  J 


(3.5) 


Substitute  t  +  hA  for  t  in  (3.2),  let  s 
index  j  to  u-h-1 .  This  yields 


-h-1  and  change  the  summation 


k 

y'  (t+hA)  =  l  a(u-h-l ,  -h-1 )y  [t  +  (u-l)AJ.  (3.6) 

u=l 


It  follows  from  (3.5)  and  (3.6)  that  y^p+1^(t)  = 

b.-^q+1^  y[t+(j-l)A],  where  b.(q+1)  =  J  b  ^a(j-u, 

J  u=l  u 


-u),  i.e. 
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(3.4)  holds  for  q+1  and  consequently  it  holds  for  q=l,2,...,k.  Let  B 
be  the  kxk  matrix  with  q^  row  (b^  ^  jbg^*.  • .  ,b^q^ ) .  If  B  is  non- 
singular  then  we  have  from  (3.4): 


y  (t)  =  l  <iy(J)(t), 
j=l  J 


(3.7) 


where  the  c.'s  are  the  corresponding  elements  of  B  .  If  the  rank  of  B 
is  less  than  k,  ther  there  exist  constants  d^ . . .  ,d^,  not  all  zero, 
such  that 

I  d  bd(q)  «  0, 

q=l  M  J 


(3.8) 


j=l,2,...,k.  Hence  it  follows  from  (3.4)  that 


l  dq  y 

q=i  H 


(q)  = 


o. 


(3.9) 


Consequently  y(t)  e  A  satisfies  either  (3.7)  or  (3.9)  and  hence  y(t)eA, 
i.e.  AcA. 


The  converse  inclusion  A^A  follows  from  Lemmas  4  and  5  in  Appendix 
A.  Indeed,  let  y(t)eA,  i.e.  let  y(t)  be  of  the  form  (A. 20)  as  in 
Lemma  4.  The  lambdas  in  (A. 20)  define  the  matrix  L  of  Lemma  4.  By 
Lemma  5  det  L  f  0.  Hence  there  exists  a  unique  solution  vector  b  of 
equation  (A^21 )  in  Lemma  4.  Consequently,  y(t)  satisfies  (A. 22)  and 
hence  y(t)eA,  i.e.  Aca.  This  completes  the  proof  of  the  proposition. 

Thus,  if  we  either  assume  that  tabular  data  can  be  approximated  by 
an  element  of  A,  or  that  the  derivative  of  the  corresponding  function 
can  be  expressed  as  a  linear  combination  of  certain  data  values,  then 
in  view  of  propositions  1  and  2  an  approximating  element  can  be  ob¬ 
tained  by  constructing  an  appropriate  autoregressive  model  such  as  (2.8) 
in  the  preceding  section.  In  other  words,  the  method  described  in  the 
preceding  section  yields  the  exact  derivative  for  a  function  in  the 
set  A=A,  or,  equivalently,  in  the  set  of  functions  that  satisfy  a 
family  of  difference  equations  of  the  type  (2.4). 


4.  DERIVATIVES  WITH  SIMPLE  EIGENVALUES 

Equations  (A. 21)  and  (A. 22)  in  Appendix  A  yield  derivatives  of 
tabular  functions  representable  by  autoregressive  models,  including 

those  with  simple  eigenvalues.  The  needed  weights  b.  in  (A. 22)  are 

J  “  ^ 
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obtained  by  solving  the  linear  system  of  equations  (A.21).  This  re¬ 
quires  complex  number  arithmetic  since  the  coefficients  in  the  equation 

(A. 21)  are  imaginary  for  imaginary  eigenvalues  A.. 

J 

A  closed  form  expression  in  terms  of  components  of  complex  eigen¬ 
values  can  be  obtained  for  the  case  of  simple  characteristic  roots.  An 
element  of  A  corresponding  to  this  case  is 


y(t)  =  l  C,XV  ,  (4.1) 

j=l  J  J 

i.e.  it  belongs  to  the  subalgebra  of  A  generated  by  trigonometric  and 
exponential  polynomials.  The  case  of  simple  eigenvalues  is  of  practical 
interest  for  the  following  reason. 

Equality  of  two  eigenvalues  defines  a  functional  relation  between 
the  a.'s,  the  coefficients  of  the  corresponding  autoregressive  model, 
i.e.  'it  defines  a  hypersurface  in  the  k-dimensional  Euclidean  space 

of  a-j^j . .  •  Thus,  the  set  of  vectors  a  in  that  define 

multiple  roots  has  Lebesgue  measure  zero.  Therefore  if  we  assume  that 

experimental  data  yields  points  of  aeX^,  all  equally  likely  in  a 
finite  region,  i.e.  with  probability  distribution  proportional  to  the 
Lebesgue  measure,  then  the  probability  of  obtaining  equal  or  nearly 
equal  eigenvalues  is  very  small,  although  synthetic  data  derived  from 
algebraic  polynomials  does  yield  multiple  roots. 

In  order  to  obtain  a  closed  form  solution  of  equation  (A.21)  we 
need  the  following  lemmas. 


(j) 


Lemma  6.  Let  Vn+^  be  the  Vandermondian  of  Xq,x-| 
be  its  minor  corresponding  to  the  element  x^. 


» ^2 , 

Then 


,xn  and 


V  ^ J )  =  s  V  (n} 
n  n-j  vn 


(4.2) 


where  S  .  is  the  symmetric  function  of  x,  ,x.,,...,x  of  order  n-i  and 
(n)  .  J  l  d  n 

Vn  is  the  Vandermondian  of  x^,X2»...,x  . 


Proof.  By  expanding  Vn+1  with  respect  to  the  first  row  (the  row  con¬ 
sisting  of  powers  of  xQ)  we  get 


21 


Divide  both  sides  by  (-l)n  V 

n 


V  ,  n  .  .  V  0) 

(-Ijn  -Q+Lj-  =  ^  (-,  )J+nx( J  — 


V  1ST  -  <£„  A0  (a)  • 

n  J  n 


(4.3) 


i  \ 

Since  V  '  '  is  itself  a  Vandermondian  we  have 


(n). 


n  j-1 


n  J-1 


v  w  -  n  n  (x.-x.).  Also,  vn+1  =  n  n  (x.-x.). 

n  j=2  i  =1  J  1  n  1  j=l  i =0  J  1 


Therefore  we  get  from  (4.3) 
n  n 

j 


I  U0  -  X.)  =  l  (-l)J  +  n  Vx 

i=i  u  J  j=o  u  v/n; 


U) 


(4.4) 


The  right  hand  side  of  this  identity  is  thus  a  polynomial  in  Xq  of 

degree  n  with  the  roots  x-j  . xn>  and  the  coefficient  of  Xq11  is 

equal  to  one.  Therefore  the  coefficient  of  XqJ  is  (-l)’1-"1  Sn_j.  Thus, 

we  have  from  (4.4)  Sn  .  =  which  is  (4.2). 

Lemma  7.  Let  S.  be  the  symmetric  function  of  A, , A9 , . . .  ,A,  of 
order  i  and  S^  p;  be  the  symmetric  function  of  .Ag,. . .  »A  -j  »X  +-| ». . 
of  order  i .  Then 


si(p)  -  jo  (-1,J  si-/pJ 


(4.5) 


Proof  by  induction  on  i.  For  i=0  (4.5)  is  an  identity  1=1.  Suppose 
(4.5)  holds  for  i  <  k.  The  function  S^+-j  can  be  expressed  as  the  sum 

of  all  the  products  of  i+1  factors  taken  from  A-|  .Ag, . .  ,A  -j  ,A  +-| , . . .  ,A^ 

plus  the  sum  of  all  the  products  of  i  factors  taken  from  the  same  set 

with  each  product  multiplied  by  A  ,  i.e.  S^+^  =  S.+-|^  +  S^^Ap  or 

S1+1(P}  =  S^+-|  -  S..^A  .  By  substituting  S^  ^  from  (4.5)  we  get 


(P)  . 


Si+1 =  Si+1  "  E  (_1)j  Si-j  X 

J 


j+1 


J  P 


Let  j=u-l .  Then 
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i+1 


(p) 


Si+1 


i+1  , 

-  I  (-1)U',S, 


U=1 


i-u+1 Ap 


i+1 

l  H)US-+1  ,/n-  This  is  the 

u=0  1+1 "u  p 


relation  (4.5)  with  i  advanced  to  i+1  and  the  proof  is  complete. 


When  all 

eigenvalues  are 

simple  equation  (A. 21) 

is  of  the 

form 

/,  s+1 

A1 

,  s+2  .  s+k' 

A-|  •  •  •  •  A-| 

\ 

bl\ 

f  r  \ 

/ log 

.  s+1 

A2 

s+2  ,  s+k 

A^  •  •  •  • A^ 

\ 

»  i 

j  j 

• 

b  ' 

j 

log  A 2  { 

1 

Ar 

A  5+1 

vAk 

.  s+2  s+k 

Ak  ’  *  "Ak 

}  1 

.  1 

/  I 

-  / 

lvbi/ 

i--- 

y°!i \/ 

(4.6) 

where  A  is  the  spacing  used  in  determining  the  A's  (i.e.,  q  times  the 
original  data  spacing).  The  solution  b  of  (4.6)  is  obtained  by  Cramer's 
rule.  This  allows  us  to  cancel  common  factors  in  the  numerator  and 
denominator  that  are  differences  of  eigenvalues  and,  hence,  to  reduce 
the  round-off  error. 


The  determinant  of  (4.6)  can  be  written  as  follows: 


det  L 

-  s  3+1  vk  =  n  a  3+1 

k  k  j=1  J 

n 

j>i 

(A. -A.). 

J  i 

14.7) 

Let 

k 

my 

=  (vv 

(4.8) 

17p 

and 

k  L  ■ 

P(x)  = 

I  a  Ak-J 
j=0  J 

(4.9) 

(ao 

=  -1 )  be  the 

characteristic  polynomial 

of  (3.1 ) .  Then 

1L 

9A 

k 

=-  n 

X=  A  j=l 

( A  -  A.).  Therefore 

r  J 

p  m 

k-1  .  . 

«V 

=  -  I  (k-j)  a.A  k_J 
j=0  J  p 

“1 

(4.10) 

or 
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(4.11) 


Q(XP)=-J1  "VqS 


q-i 


By  writing  n  (X.-X.)  =  n  (Xn  -  X.)  n  (X.  -  X.)(-l)k“p  and  by 
j>i  3  i=l  p  1  j>i  J  1 

i^P  i  »j7p 

observing  that  the  last  product  is  a  Vandermondian  of  X,  ,XOJ...,X  , 

^  /_\  id  P“ I 

Vl’‘"’V  say  V^_^p  ,  we  can  write  (4.7)  as  follows: 


det  L  =  (-1  )k_p  Sks+1  Q ( X  ) 


(P) 


(4.12) 


Continuing  the  application  of  Cramer's  rule,  let  s  be  the  matrix 

obtained  from  L  by  replacing  its  (j-s)-th  column  by  the  right  hand  side 
of  (4.6).  By  expanding  det  B.  with  respect  to  this  column  we  get: 

(4.13) 


det  B . 
J-s 


.  (-l)j 


-s 


S+l  _s_l  ( j-S-1) 


A' 


J,  <-up  i-gVrx-X,  <’>• 


where 


k-1 


(j-s-1 ) 


(p)  = 


1 

A1  • 

C\J 

1 

CO 

I 

r<~ 

>  J-s 

A1  • 

X  k“1 

..x1 

1 

Vi- 

3 j-s-2 

••Ap-1 

J-s 
AP-1  “ 

.  x^1 

P 

1 

Vi  ■ 

y J'-s-2 

••Vi 

J-s 

Vi 

J-i 

Vi 

1 

\  ■ 

*j-s-2 
*  ’  Ak 

J-s 

Ak 

•  i 

1 

7T  7T 

1  1 

1 

(4.14) 


i.e.  vVs"Vp)  is  a  corresponding  minor  of  det  B.  divided  by 
k-i  j-s 

cs+l,-s-l 
bk  Ap 


Determinant  (4.14)  is  a  minor  of  a  Vandermondian  in  X^,X-j,...,X  p 
Xp+p...,Xk  corresponding  to  the  (j-s-1  )-th  element  of  the  first  row. 
Thus,  by  Lemma  6  above  we  have: 
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V  (J-S-1)/  x  =  s(p)  u(p) 
vk-l  lp;  Vj+s  vk-l  ’ 


(4.15) 


where  sfp}, 
k-j+s 

Ap-1  ,AP+1 


is  the  symmetric  function  of  order  k-j+s  in  X,  ,X0,...» 

~(p)  1  4 

,Xk  and  Vk_- j  is  the  Vandermondian  of  the  same  arguments. 


By  substituting  (4.5)  in  (4.15)  we  get 


v( j-s-1 ) 
vk-l 


(P) 


(-1)V  S 


k-j+s-v  p 


(4.16) 


Since  S.  is  a  symmetric  function  of  the  roots  of  (4.9)  we  have 


(-1) 


k- j+s-v-1 


k-j+s-v  ~  v-,/  “k-j+s-v  ' 

By  substituting  (4.15)  -  (4.17)  in  (4.13)  we  get: 
k 


(4.17) 


det  B 


1 


Ar  p=l 


\  H)PlogVsk+1  XpS"1Vk-l 


k-1  .  v 

k-j+s-vAp 


l  (-1 )"  ‘  a,.  ..X 


k-j+s 

*  l 

v=0 

In  view  of  (4.11),  (4.12),  and  (4.13)  Cramer's  rule  yields 

k-j 

v=0  ak'J'“v  "P 


(4.18) 


i,  X  s  logrX_  y  a.  .  Xv 

b  _1_  5  P  py.Q  k-J‘v  t 

3  Ar  pSl  *  .v 

l  vXn  a. 


(4.19) 


v=l 


p  k-v 


Formula  (4.19)  provides  a  closed  form  expression  for  the  coefficients  in 
(3.2)  in  terms  of  the  coefficients  a^  of  the  autoregressive  model  (2.9) 

and  in  terms  of  the  eigenvalues  of  the  corresponding  difference  equation 
Since  some  of  the  eigenvalues  may  be  imaginary,  computation  of  the  b.'s 

by  (4.19)  still  may  require  complex  number  arithemetic.  This  can  be 
avoided  by  combining  complex  conjugate  terms  in  (4.19). 

^  ^ 

Let  ( Xj ,X| ) ,  (X,,X.  (XU,X)  be  the  pairs  of  complex  conjugate 

roots  and  ^u+l 5 * •  •  rea^  roots  of  the  characteristic  polynomial. 
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Write  (4.19)  as  follows: 


bj 


1 


u 

F  I 

A  p=l 


PS  l091"  AP  in  ak-j-v  XP 


k-j 


v=0 


V  v 

V  ak-v  Xp 


+  T 


1 


A  p=2u+l 


*  APS  ^  Ap  Y  Vj-»  i 


v=0 


l  v  ak-v  xp 
/  =  1 


where  T*  is  obtained  from  the  first  term  in  the  brackets  by  replacing 
A.  by  A  .  By  adding  the  terms  in  the  brackets  we  get 

r  r 


( 


,  r  n  k;j  v\  k  *v  /  *-S  r  *  k;j 

AP  105  AP  Jo  aM- v  Ap  j  i,dvxp  +(*p  '"S  \  "k 


*v\  k  v 

•  »A  )  7  d,.A 
-J-v  P  /vSj  v  P 


(4.20) 


where  d  =  v  a, 

v  k-v 


The  denominator  Dp  of  the  last  expression  can  be  transformed  as  follows: 


DP  ■  j,  4<VP*)V  +  vIn  <W  (x/xj"  ♦  X"pxf, 


By  writing  A  =  A  +iB  ,  a  -  log  | A  | ,  and  8  =  arctan  B  /A  we  get 

rrrr  r  r  rr 

k  k  (v+n)a  r  -i 

D  =  T  f  e  p  cos  (v-n)8  d  d 
P  V=1  n-1  L  pJ  vn 


Similarly,  the  numerator  Np  of  (4.20)  can  be  transformed  into 
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k- j  k 

NP  =  2  l  F  ak_ i_ V  dn  exP 

p  v=o  n=1  K  J  n 


(v+n-s)a  +  r  log  Vet  2+3  ^ 
P  P  P 


COS 


(v-n-s)3p  +  r  arcsin^Bpy^c*p2+3p2  j 


Neither  N  nor  D  involves  any  imaginary  elements. 

r  r 


Consequently,  we  get: 


k-j 


lL  y  V  1_  *  V  '°9  *pJ«Wp 

J  Ar  p=l  up  Ar  p=2u+l 


v=0 

~k  7 

r, v  ak-vxP 

v=i  K 


a  closed  form  solution  for  the  coefficients  in  (3.2). 


(4.21) 


5.  ERROR  BOUNDS 

In  this  section  we  develop  expressions  which  yield  approximate 

bounds  on  the  errors  in  derivatives  calculated  by  the  method  of  Section 

4. 


ket  x^.  =  x(t),  t=0,A,2A, . . . ,nqA  be  tabular  data  of  a  function 
y(t)  eA  with  the  measuring  error  e^.,  i.e. 


x(t)  =  y(t)  +  et.  (5.1) 

We  assume  that  is  a  white  noise  with  zero  mean  and  variance  a2.  As 

shown  on  the  preceding  pages  the  r-th  derivative  of  y(t)  is  expressible 
in  the  form 

(rl  t+s 

y  ;U)  =  I  b.  (r)  y(t+jqA),  (5.2) 

j=s+l  J 

and  its  approximation  x^  (t)  is 


X 


(r) 


(t) 


k+s  ^ 

l  b.  (r)  x ( t+jqA) , 
j=s+l  J  6 


(5.3) 
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where  the  b.  's  are  expressed  in  terms  of  the  coefficients  in  an  auto- 
J  ~  S 

regressive  model 

k 

x(t)  =  l  a.  x(t-jqA)  (5.4) 

j=l  J 

corresponding  to  the  difference  equation 
k 

y(t)  =  l  a.  y(t-jqA)  (5.5) 

j=l  J 

satisfied  by  the  function  y(t).  Let 

<S(r)  (t)  =  x^(t)  -  y^(t)  (5.6) 

/\ 

denote  the  error  of  the  approximation.  Since  the  coefficients  b. 

J  “  S 

are  functions  of  tabular  data,  we  can  use  the  relations  (5.1 )-(5.6)  to 

express  6V  ' (t)  as  a  function  of  x(t) ,t=0,A, . . . ,  and,  thus,  obtain  the 

variance  of  <S^(t)  expressed  in  terms  of  statistics  of  This  yields 

a  measure  of  the  accuracy  of  numerical  derivatives  obtained  by  the 
method  described  in  the  preceding  sections. 

In  the  present  section  we  derive  an  approximation  to  the  bound  of 
(r) 

the  variance  of  6V  '(t) .  To  this  end  let 

ej-s  =  Vs(r)  "  bj-s(r)  ’  (5*7) 

i.e.  3-  ~  is  the  error  in  the  estimate  of  b.  (r).  In  order  to  simplify 

3 J “S  ^ 

subsequent  algebraic  expressions  we  assume  that  the  estimates  b._  (r) 

J  s 

are  unbiased,  i.e. 

E  (3.  J  =  0.  (5.8) 

By  substituting  (5.2)  and  (5.3)  in  (5.6)  we  get 
<5(r)(t)  =  At  +  Bt, 

where 

k+s 

At  =  J  Jj  bj-s  et+jqA 


(5.9) 


(5.10) 
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and 


k+  s 

Bt  l  sx(t+J'qA)- 

j=s+l  J 


(5.11) 


It  follows  from  (5.8)  -  (5.11)  that  E(S^ (t) )=0.  Therefore  by  squaring 
both  sides  of  (5.9)  and  taking  expected  values  we  get 

oZ{r)  =  E(A^)  +  E(B^)  +  2E(A^.Bt) .  (5.12) 

2 

where  a  (r)  is  the  variance  of  the  error  in  the  r-th  derivative. 

Since  E(B^At)  <  \Ze(A^)  E(B^)  we  get  from  (5.12): 


o  (r ) 


<(^  +/e(B^)  ) 


(5.13) 


Thus,  E(A^.)  and  E(B^)  define  a  bound  on  a  (r). 
It  follows  from  (5.10)  that 
o  9  k+s  9 

E(At)=a  l  b  s(r)' 

j=s+l  J"S 


(5.14) 


We  replace  the  values  b.  (r)  in  (5.14)  by  their  estimates  and  obtain 

J  « 

an  approximation: 


O  o  k+S  An 

E(A  )  =  o  l  b2.  (r). 
j=s+l  J"s 


(5.15) 


Let  Aj  and  X . ,  j=l,2,...,k  be  the  eigenvalues  of  the  difference  equa¬ 
tions  (5.5)  and  (5.4),  respectively.  Let 


and 


a.  =  a  .  -  a . 
J  J  J 


nJ  =  XJ  "  V 


(5.16) 


(5.17) 


We  assume  that  all  the  eigenvalues  of  (5.4)  and  (5.5)  are  of  multi¬ 
plicity  one.  Then  the  coefficients  b.  (r)  and  b.  (r)  are  given  by 

J- S  J “S 

.  /\ 

(4.19)  in  terms  of  A.'s  and  a.'s  and  A.'s  and  a.'s,  respecti vely .  An 

J  J  J  J 
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expression  for  Bt  in  terms  of  the  ou's  can  be  obtained  from  (4.19), 

(5.7),  (5.16),  and  (5.17).  This  is  a  non-linear  al gebraic . function. 
We  approximate  this  function  by  the  following  linear  relation 


k+s 

k 

dB, 

l  x (t+jqA) 

l 

J- 

dry 

s 

ai 

j=s+l 

i=l 

U  Lvj 

k 

k+s 

dB, 

Bt  =  l 

1  i=l 

l 

=s+l 

_ 1 

d 

or 


x (t+jqA) , 


where 


(5.18) 


dB.  e 
J~s 

da. 


dB . 
P3-s 

da. 


k  3B  • 

3A„ 

=  I  J_S 

ai =0  P=1  8Xp 

*  _ E. 

/N 

A  =A 

P  P 

3a^. 

a 


i=0 


33, 


J-s 


3a 


ai  =0 


and 


3B  • 

Jzl 


3a  1- 


a. 


i  =0 


3B- 

J-s 

3a 


a .  =a  • 

l  i 


(5.19) 


In  order  to  simplify  further  derivation  of  an  expression  for  E(B^  )  we 
introduce  the  following  notation: 


Np  -  A'2  logr  Xp, 

k-j+s  u 

Npj  =  UJ  ak- j+s-u  y 


‘  i,  u  xp  vu* 


and 


N( j,p)  =  N  N  . 
Vk>  Vl  P  PJ 


ffx  \  \  a  a  a  1  =  — - ^ *  P - 

T\A1  »a2»  •  •  •  »Au»di  »a9».  •  • »“!/'  n 


k  *  1  2 : 


(5.20) 

(5.21) 

(5.22) 

(5.23) 

(5.24) 
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Then  it  follows  from  (4.19)  and  (5.7)  that  3.  =  f(A-, 

J  s 

-  f ( X-j ,. . .  ,Xk>a1 . .  ,ak) .  Therefore, 


3ei-s 

3A_ 


D  dN  (j  ,p)  ^/ .  v  — £ 

Up  3A  ^  v  J  , P /  ax 

— _ p _ e_ 


A=An 

P  P 


'  \»al  ,ak  ^ 


(5.25) 


By  (5.23) 


AN(j,p)  a  ?0sL  N  +!!k  N 

9Xp  PJ  9\>  v 


(5.26) 


By  (5.20) 


9N  sN_  rN_ 

_ R  = _ E.  +  P  , 

9A„  A  A  1 ogA„  ’ 

P  P  P  P 


(5.27) 


By  (5.21) 


9N  •  k"J+s 

'  i  u  ak-j+s-u 


u-1 


3A 


p  u=0 


(5.28) 


By  (5.22) 


3D. 


~  l  u  A 


2  .u-1 


p  u=l 

Similarly  we  get 


3A_  A,  VP  “k-u 


awtJ.p)  n  .  !V(i  . 

k  3a,  up  3a7lu,p' 

1  2 

p=l  D„ 
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j-s 


9ai 


a .  =a  • 

l  l 


(5.29) 


(5.30) 


where 


=  "n 

9a  n-  P  P 


(5.31) 


and 


9D 

3i^  =  (k_i)Ap 


•x,  k-i 


(5.32) 


k  k-u 

By  definition  of  Ap  we  have  Ap  =  l  auAp  .  We  differentiate  this 
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expression  with  respect  to  a.  and  solve  the  resulting  equation  for 

3V3ar 


3d  . 

1 


k 


(k-u)a  A 
v  up 


k-u-lj 


(5.33) 


The  formulas  (5.19)  -  (5.33)  define  the  coefficients  d£.  /da.  in  (5.18) 

J  “S  1 

and,  thus,  the  quantities  in  the  following: 


k+s  dB . 

ci  =  l  -foT1  x ( t+jqA) .  (5.34) 

1  j=s+l  dcti 


k  k  k 

We  have  by  (5.18)  B  =  £  C.  a,.  Hence  B  2  =  l  l  a  a  C  C  and 

1  i=l  1  1  r  v=l  u=l 


E(Bt2)  =  Et[F.(Bt2|t)] 


1  CuCvE(Vv|t> 

V=1  U=1 


(5.35) 


where  the  conditional  expectation  is  over  the  distribution  density  func¬ 
tion  of  the  a. 's.  Since  E(auav|t)  does  not  depend  on  t  we  get 


E(Bt2)  ■  Z  I  E(aucv).  Et(CuCv)  ■  \  Z  E(Vy)  •  J  I  Cu<t)Cv(t) 
V  =  1  U=1  v  =  l  u  =  l  t 


1 

N 


k  k 

l  l  l 

t  v=l  u=l 


E(“u%)Cu(t)Cv(t)- 


^  2 

The  quantity  No  in  (2.17)  is  small  as  compared  to  the  diagonal 
p  ^ 

elements  of  the  matrix  M.  Therefore  the  estimates  a.  in  (5.4)  differ 

J 

only  little  from  the  least  squares  estimates  of  the  a.'s,  and  hence 

J  2  ~  2  -1 

the  covariance  matrix  (E(aua^))  can  be  approximated  by  a  (M-Ncz  I)  = 
cr2M  .  By  writing  m  for  an  element  of  M  obtained  by  estimating  the 
regression  coefficients  (a.'s)  we  have  the  following  approximation: 

J 

E(B2)  =  {£  l  Z  I  C  (t)Cv(t)m  (5.36) 

r  IN  t  u=l  v=l  u  uv 
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aj's»  formula  (5.36)  yields  an  estimate  of  E(Bt2).  This  together  with 

(5.13)  and  (5.15)  defines  an  upper  bound  of  the  variance  of  the  error 
in  the  estimate  of  the  r-th  derivative.  However,  this  is  not  an  absolute 
bound  since  its  estimate  depends  on  the  assumption  that  (4.19)  yields 
unbiased  estimates  of  the  b j ' s ,  on  the  accuracy  of  linearized  expressions 

of  bj  s  in  terms  of  a ^ 1 s  and  X j 1 s ,  and  on  the  effect  of  replacing  a.'s 

and  Aj  s  by  their  estimates,  a.'s  and  A.'s,  respectively. 

A  few  examples  listed  in  Appendix  B  show  that  the  formula  (5.13) 

yields  reasonable  bounds  of.a2(r)  in  almost  all  the  cases  tested  and 
hence  it  can  be  used  as  an  indicator  of  accuracy  for  the  numerical  esti¬ 
mates  of  derivatives  obtained  by  the  method  of  this  report. 

Further  applications  and  examples  will  be  presented  in  a  companion 
report  now  in  preparation.  Additional  sets  of  synthetic  data  as  well  as 
samples  of  field  test  data  will  be  considered  therein. 
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appendix  a 


Lemmas  1,  2,  3  below  simplify  proofs  in  Section  2,  while  Lemmas 
4,  5  are  found  useful  in  Section  3. 


Leiyma  b  Suppose  y(t)  satisfies  the  conditions  (a),  (b),  and  (c) 
of  Section  2.  Let  S^.(A)  be  the  symmetric  function  of  the  roots  of  (2.4) 

of  order  j  and  M^a)  =  Sj(A)  -  Qj  .  Let 
k  .  . 

N(A)  =  l  (-l)J_i  M.(A)y(t-jA) .  (A.l) 

j=l  J 


Then 


dqN(A) 

dAq 


A=0 


=  0 


for  q=0, 1 , . . . ,k-l . 


(A. 2) 


Proo_f.  With  f(t)  replaced  by  y(t)  equation  (2.4)  can  be  written 
as  follows  by  expanding  PA,  replacing  powers  of  B  ,  and  solving  for 

y(t) : 

k 


y(t)  =  l  (-l)0’"1  S  . ( A )  y(t-jA)  . 
j=l  J 

By  substituting  (A. 3)  into  the  identity 
0-B.)ky(t)  =  y  +  |  (-1)J  k 

j=l 


(A. 3) 


t(t-jA)  we  obtain 


(l-BA)ky(t) 


I  (-l)J'1M,(A)y(t-jA)=  N(A). 
j=l  J 


It  follows  from  (A. 4)  that  the  k-th  derivative  of  y(t)  is 


(A. 4) 


y(k)(t)  =  lim  .  (A. 5) 

A^O  AK 

Now  we  prove  (A. 2)  by  induction  on  q.  Indeed,  it  follows  from  the  def¬ 
inition  of  Mj (A)  that  Mj (0)=0  since  each  term  of  S^O)  =  1.  Hence  by 

(A.l)  N(0)=0,  i.e.  (A. 2)  holds  for  q=0.  Suppose  it  holds  for  q=n<k-l . 
Then  by  applying  L’Hospital's  rule  to  (A. 5)  n+1  times  we  get 
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,<k)(t)  .  11B  dn+IN(A)/dAn+1 


A-*0  k(k-l) . . .  (k-n)A 

,00 


k-n-1 


(A. 6) 


By  hypothesis  (b)  y'  ;(t)  is  bounded.  Therefore  (A. 6)  implies  that 


dn+1N(A)/dan+1 


n  =  0,  i.e.  (A. 2)  holds  for  q  =  n+1  <  k.  Q.E.D. 


A=0 


For 


j-1 


nj  and  X.  in  (2.2)  let  =  J  (n.+l)  and  vm  +1  =  vm<+2  =  . 

i  -1  J  J 


v 


=  A..  Let  T,  .  be  the  set  of  all  the  combinations  of  the  integers 
j+1  J  KsJ 

l,2,...,k  taken  j  at  a  time.  Let  M.  =  M . ( A) . 

J  J 

Lemma  2.  Let 


L  (i  i  ...,i  )  =  l  log  v 
K  1  J  u=l  \ 


Then 


(A. 7) 


Z  (-DJ  JP  I  Lj-P  (i  1  )  »  0 

J=1  TkJ 

for  p  <  q  <  k. 

Proof.  Differentiation  of  (A.l)  yields 


But 


M. 

J 


(q-p)  . 


(A. 8) 


Z  (-1)-3'1  [M.y(q)(t-jA) 
dAq  j=l  L  J 

♦  Y  (jj)  Mj(c,-p)  y(p)(t-jA)  (-j)pJ  . 


(A. 9) 


l  v.  A  V.  A...v.  ALq-p  (iri2,...,i .), 

Tk#J  ’l  n2  nj  k  1  2  J 


36 


Thus , 


dqN(A) 

dAq 


A=0 


-  I  (-1)3'1  V  y(p)(t)  (-j)P 

j=l  p=0 


k,  j 


q-p 


( 1 1 5 1"  2 5  • 


’V 


By  interchanging  the  order  of  summation  we  get 


dqN(A) 

dAq 


-V 

A=0  p: 


E1  yCp)(t)  f 2]  i  (-i)J-i-PjP 

=0  \pj  j-1 


•  £  Lq_p  ( i  i  i  \ 

T  Lk 

k,j 


(A. 10) 


It  follows  from  (A. 2)  and  (A. 10)  that 


V  y(p,(t)  (M  I  (-I)J-I-P  jP  j  L^d  i . 1.)=0 

p=0  \  P I  1=1  T  Kid  j 

k,j 

(A. 11) 


for  q  <  k. 


Let  pQ  <  q-1  be  the  largest  value  of  p  with  non-vanishing  coeffi¬ 
cient  of  in  (A. 11),  i.e.  (A. 11)  can  be  written  in  the  form 


(Pq)  Pr  1  ~  (p) 

y  +  L  C  y  =  0,  an  ordinary  differential  equation  of  order  p„ 

p=0  p  K0 

m  "j  .  V 

whose  solution  y(t)  is  thus  of  the  form  y(t)  =  l  l  C-.t1  e 

j=l  i =0  J1 
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with  pn  =  l  (n.+l).  But  it  then  follows  that  y(t)  satisfies  the 

U  j=l  J 

m  u .  n  .+1 

difference  equation  P(B)y(t)=0  where  P(A)  =  £  (A-e  J  ,  i.e.  a 

j=l 

difference  equation  of  order  Pq  <  q-1  <  k,  contrary  to  the  assumption 
(c).  Therefore  in  (A. 11)  we  must  have 


I  (-DJJP  l  Lj-p  (i  — ,ii)=o, 

J  k,  j 


(A. 12) 


for  p  <  q  <  k,  which  is  (A. 8)  .  Q.E.D. 

Next  we  prove  the  following  lemma,  making  use  of  the  previous  one. 
Lemma  3.  Let  •  •  •  *i j)  be  as  in  (A. 7).  Then 


k,u 


(■!)" 

u!  (k-u) ! 


i  (-nJjk-u  i 

Tk,j 


•V'v 


(A. 13) 


for  u=l,2,...,k,  where  Sk  is  the  symmetric  function  of  £.j=log  , 
i=l  ,2 , . . . ,k,  of  order  u . 

Proof.  This  lemma  is  proved  by  induction  on  k.  When  k=l  the 
right  hand  side  of  (A. 13)  is  -1*(-1)*£^  which  is  j.  Suppose  that 

(A. 13)  holds  for  k-1.  For  any  positive  integer  u,  S.  is  a  linear 

K  y  U 

function  in  each  argument  and  it  can  be  written  as  follows: 

^k,u^k^  =  ^k  ^k-1  ,u-l  +  Vl.u 


for  u=l,2,...,k  (we  assume  that  S 


k-1  ,k 


(A. 14) 

=  0  and  S\  g=l  for  k  >  1).  Since 


V  in  (A. 13)  is  a  symmetric  function  of  the  E. 's,  we  need  only  to  show 
that  V  is  representable  in  the  form  (A.  14),  i.e.  that 


?k=0  =  Sk-l,u 


(A. 15) 
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and 


9V_  =  ? 

Kk  Vl.u-l  • 


(A. 16) 


We  write  (A. 13)  as  follows: 


V(?k}  =  uT(k-u) ! 


i-1.1  L*-i<0  + 


+  (-l)k  kk-u  Lk_x  (1 ,2, . . .  ,k-l ) 


k-1  .  . 

+  I  (-DJ  Jk'U 

j=2 


3k-U  [I  L^_ 
L  k-1 ,  j 


1  ( 1  1  ^ ) 


Therefore, 


+  I  ^  Lk- + 

k-1  ,j-l 


V(0)  = 


(-1)K 

u [ (k-u) 


.  (A. 17) 


-  1  LjJ.jd)  +  (-l)kkk-Vk_l(l,2 . k-1) 

Tk-1 ,1 


+  1  (-1)J  jk"u  I 
J=2  Tk-l,j 

k-1  .  . 

+  l  (-1)J  jk"u  I  Lk-l(1l,i2”  ’  ‘  »i  j_i ) 
3=2  Tk-1 , j-1 


39 


By  changing  the  summation  index  in  the  last  sum  and  by  combining  the 
terms  we  get 


V(0)  = 


5T{^r|(2k-u-l)  l  LjJ.jd) 

'  k-1 , 1 


+ X  (-dj  [jk'u  -  <j+i>k'u]  •  i 

3  k-l ,  j 

(1>2 . M)J. 

binomial  theorem  to  the  expressions  in  the  brackets  and  get 

l  Lk_i  (■>  i » •  •  •  »1j) 


-1J) 


We  apply  the 


(  i  (  k”l  •  ,  , 

v(°)  ■  -sTriUr  j  Ij  (-DJ(K-u)Jk-u-1 


Tk-1 ,  j  f 


j=l 


i  =2 


k-1  .  k-u  /k-u  ,  , 

M.H)'  lltji1'1-'  Z  L“.1o1,...ij)  . 


T, 


k-lj 


By  the  inductive  assumption  the  first  term  in  the  braces  multiplied  by 
the  factor  outside  the  braces  is  ^  u<  We  interchange  the  order  of 


summation  in  the  second  term  and  get: 

t k-1  k-u  /k-u 

I 

i  =2 


_  /  -i  \  k- 1  k-u  / k-i 

V(0)  -  S.  .  +  Kfr — rr  J  \  i 

k-l,u  ul(k-u)!  . L0  x 


k-1 
j 


Zt  (-DJ  Jk'u'1  Z 


k-1 ,  j 


By  (A. 8)  the  expression  in  the  brackets  is  zero,  i.e.  V(0)  =  S,  , 

k-l,u 

Thus,  (A. 15)  is  proved. 
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In  order  to  verify  (A. 16)  we  differentiate  (A. 17)  with  respect  to 
£k  and  get: 


9V  (-1)K 
Kx,  u! (k-u) 


r  j-u^"1  +  u(-l)k  kk'u  |^Ltc_1(l ,2 . k-l)+£; 


u-1 


u  [Lk-l^r*  ’•  ’Vl^k 


k 

u-1 


+  X  ^  ^  T  I  "Mc-l'T 

J  k-l,j-l  L 

(A. 18) 

By  applying  the  binomial  theorem  to  the  expressions  in  the  brackets 
and  by  combining  the  terms  containing  the  same  powers  of  we  get: 


9V  _  (-1)K 

9L  (u-1)  1  (k-u)  1 


•V  (-i)3+1(j+i)k-u  i  l“;{  dj . i  ) 

Tk-i,j 


.1  +  i  (-i)J  /-u/k'1' 
j=2 


.u-1 


d  I  (-Dj  jk"u  I 


j=2 


T 


k-1 ,  j-1 


Lu-i-i  (1 
Lk-1  ul’ 


(A. 19) 


By  applying  the  binomial  theorem  again  and  by  interchanging  the  order 
of  summation  the  coefficient  of  the  zero  power  of  in  the  braces  can 

be  written  as  follows: 


k-u 


-  I 

i  =0 


k-1 


(-l)j 


.  k-u- i 

J 


l 

Tk-l,j 


By  (A. 8)  all  the  terms  of  the  outer  summation  except  the  term  corre¬ 
sponding  to  i=0  vanish.  Hence  the  coefficient  of  the  zero  power  of 
5k  in  (A. 18)  is 


(u-1 ) (k-u) 1 


k-1 

l 

j=l 


(_l)j  jk-u 


I 

Tk-1 ,  j 
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By  the  inductive  assumption  this  is  -j  i .  Thus,  in  order  to  com¬ 
plete  the  proof  of  (A. 16)  it  remains  to  show  that  the  coefficients  of 
every  non-zero  power  of  in  (A. 19)  vanish.  Since  the  highest  power 

of  in  (A. 19)  is  u-1  we  need  to  consider  only  the  case  of  u  >1. 

The  coefficient  of  can  be  written  as  a  constant  times  the  form: 


k  .  .  /  k-1 

^  <-i)j  H-i 


k-l 

l  (-1)’  (1+1) 
1=0 


k-u 


-(1-B)k_1  X 


k-u 


x=k  ‘ 


u-1 


Since  u  >  1  it  follows  that  the  coefficient  of  is  zero.  The  co¬ 
efficients  of  i=l  ,2, . . . ,u-2  in  (A. 19)  vanish  because  of  (A. 8).  This 
completes  the  proof  of  (A. 13). 


-1 


We  use  the  following  notation  in  the  next  lemma:  II  (r-v)=l, 

v=0 

AJp  ■  (Wl>P  *q+I.  (S+2>P  X=+2,.....(s+k)P  X^), 


P-1 


=  — — ■  (logr-':)A  n  (r-v)  )  for  p=0,l  ,2, . . .  ,n  ,  q=l  ,2, . . .  ,m, 
Ar  \  q  v=0  /  4 


and  L  -  (A10»An>  *  *  *  »Alni ,A20»A21  ’  *  *  *  ,A2n2’  *  "VnJ  ’ 


Lemma  4.  Let  y(t)eA,  i.e. 


m  nj  i  t 

y(t)  ■  X  l  _  cj,  t1  x* 


j=l  1=0 


(A. 20) 


with  k  =  l  (n.+l).  If  there  exist  b  =  (b,  ,b9, . . . ,b.  ) 
j=l  J  1  £  K 

Lb  =  A 

then  for  every  integer  s  satisfying  -t-A  <  sA  <  T-kA 


such  that 


(A. 21) 
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(A. 22) 


( r )  s+K 

y  =  I  b-:_s  y(t+jA). 

j=s+l  J  s 

Proof.  Substitute  (A. 20)  on  both  sides  of  (A. 22)  and  equate  the 

coefficients  of  the  terms  t1'^  on  the  opposite  sides  of  the  result. 

Then  (A. 21)  is  obtained.  The  lemma  is  proved  by  reversing  these  steps 
as  follows.  K 


The  expression  (A. 21)  is  a  system  of  equations 
T  1  p~  1 

Aqp  b  =  “F  lo9r"P  \  P=0,1 . nn,  q=l  ,2, . . . , 


v=0 


m  or 


.1  <=+j)P  bJ  =  V  1o9r'P 

J=1  4  J  Ar 


P-1 

A  n  (r-v) 
q  v=0 


(A. 23) 


For  each  q-l,2,...,m  and  z=l,2,...,n  multiply  the  equations  of  (A. 23) 
corresponding  to  p=u-z  (u=z,z+l , . . . ,n  )  by  c  ^  and  sum  on  u.  This 
yields:  q  '  ' 


u- z  C()u(")  j?!  (s+J)U 


2  C1  b . 
q  j 


Ar  1  Cqu(z)l09r"U+Z  Aq  .L  (r_v)-  (A. 24) 


u-z-1 

n 

v=0 


u-z-1 


write  (A. 24)  in  the  form 
n 


u=z  Cqu  (^!i-2 )  Jj  (s+j)U 


-2  X^  b. 

q  j 


-  V  Iq  C„„  logr~u-z 


A  u=z 


qu 


u-z-1 

n 

v=0 


we  can 


^u-z)v"o  (u"v)- 


(A.-25) 
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Now  we  multiply  (A. 25)  by  tz  and  sum  on  z  from  0  to  n^: 


r  ^  i  Ju-z)  i  (s+j) 


-q  .Z  z 


z=0 


u=n 


^u“z  x!+j  b. 


3=1 


.  nq  nq 

=  7  I  I 

A  z=0  u=z 


tz  logr_u“z  X, 


q  \u-z 


r  \  u-z-1 

n  (u-v). 


v=0 


We  interchange  the  order  of  summation  on  both  sides  of  this  equality: 

n. 


k  q  ,  •  'q 

i  yq+J  I  c 

j=l  J  q  u=0 
n_ 


u  /  u 


qu 


l 


z=0 


t  (s+j)1 


=  i-  I 

V'  L 


Ar  u=0  qU  z=0 


I  tz  log 


r-u-z  x  /  r 


u-z-1 

n  n  (u-v). 

q  v=0 


Multiplication  of  both  sides  of  this  equality  by  X^j  yields: 

l  b,  f  C  X*+S+J  (t+s+j)u  =  f  c  (tV) 


j=l  J  u=0  qu  q 


u=0  'qu  dtr 


By  summing  both  sides  of  this  equality  on  q  from  1  to  m  we  get: 
m 


n 

m  q 


Z  b  I  I  c 

j=l  J  q=l  u=0  qU  q 


m  q  .r  „  + 

=  T  7  (c  tuXt). 

q=l  u=0  dtr  qu  q 


k 

In  view  of  (A. 20)  this  is  jT  b. 

3=1  J 


y(t+s+j) 


(t)  which  is  equivalent 


to  (A. 22). 

Lemma  5.  Denote  nil  =  11  jl,  n:=  X  j,  and  P„  =  II  x(s+^nj+^ 

3=1  3=1  "  3=1  J 

Let  the  matrix  L  be  as  in  Lemma  4.  Then 
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m  n  :  m  j-1 

J 


,  .  (n.+l)(n.+l) 

detL=p  Hn!!  n  A.J  n  n  (A.-A.)  1  J  .  (A. 26) 

1=1  j=l  J  3=2  1=1  J  1 

^22f*  Let  Ti  =  pi  9^7  and  TiP(Xj.)f(yi )  =  T.pf(y. ) 


u.=x. 
t  j 


i.e.  the  result  of  the  operator  T^p(Xj)  is  obtained  by  applying  the 

operator  T..  p  times  and  then  substituting  X.  for  p. .  In  particular, 

J  1 

Ti0)  UJffu.*)  =  f  ( A . ) .  Let  m.  =  7  (n.+l)  and 

J  J  i=1  i 


m . 


/  v  p  ”'j  1-m.  ,-l 

t(p)  =  n  n  t,  J_i  (x.)  . 

J 


n  t. 

j=l  i*m.  ,+1  1 

J  A 


(A. 27) 


Let 


* 

L  = 


s+1  s+2 


•  -Pi 


s+k 


s+1  s+2 


s+k 


Then  T(m)L  =  det  L.  This  identity  is  used  to  prove  the  lemma  by  in¬ 
duction  on  nj3  j=l ,2, . .  ,m.  Indeed,  if  ryO,  j-1 ,2 . m=k  then  by 

(A. 27)  the  operator  T(m)  =  T(k)  simply  replaces  p.  by  X.,  j=l,2,...,k, 

i.e.  the  left  hand  side  of  (A. 26)  is  a  Vandermondian  multiplied  by 


s+1 


k  j-1 


nXj  and  the  right  hand  side  reduces  to  II  XT  II  II  (X.-X.), 
J"1  j=l  J  j=2  i=l  v  J  1  ’ 

i.e.  (A. 26)  holds  for  this  choice  of  n.'s.  By  induction  assume  that 
(A. 26)  holds  fornl9n2 . np-1  with. 


P+1  p+2 


nn4-9  ~  •  •  •  =  n„  -  0 


m 


(A. 28) 


and  np  replaced  by  v  such  that  0  <  v  <  n  Under  these  assumptions 
(A. 27)  Implies  that  P 


Tfm)  =  T*(j,p,v)  T(p,v) , 


(A. 29) 
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for 


where  T  (j,p,v)  is  an  operator  that  replaces  y.  by  A.  n 

J  J-m  i-v+p-l 


j-m  j+v+l , . . . ,k  and 


raP-itvtl  T  '-"'p-r1 

T(p.v)  =  n  Ti 

i=vl+1 


T(p-l), 


(A. 30) 


with  T(p-l)  defined  by  (A. 27). 
Let 


"k 

P  (p) 


P-1  P-1 

n  n. : :  n 
i  =1  1  i  =1 


n.  • 

A.  P  , 

i  p-1 


n 

p-l>j>i 


(yv 


(ni+l)(nj.+l ) 

(A. 31) 


By  induction  (A. 26)  holds  for  the  parameters  in  (A. 28).  Hence  it 
follows: 


or 


T*(j,p,q)  T(p,q)  L*  =  P*(p)q  1  i  A.  q:A  (<l+1Hs+l) 

P  P 

•  n  as.+1  n  (x  -A.)(q+1)(ni+1)  n  (a.  -  a  )q+1 

j=p+l  J  i<p  p  j>p  J  P 

n,+l 

•  n  (A.-A.)  n  (A.-A.) 

j>P  J  j>i >p  J 

i<p-l 


T  (j,p,q)  T(p,q)  L*  =  P*(p)q“A  q:A 

r  r 


,,  q:,  (q+1  )  (s+1 ) 


•  II  /,  .  \(q+l)(n.+l)  *,  .  .  i 

(W  1  T  (J’P^)  11  (y^-Ai ) 

i<P  v  j>m„  ,+q+l  J 


n.+l 


j>m  i+q+i 
i<p-l 


n  (y-yJ 
J>i>mp_i+q+1  J’>mP-i 


n  yf1  n  (y  -A  )q+1 
+q+l  J  j>p  J  P 


(A. 32) 
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Since  T  (j,p,q)  has  an  inverse  it  follows  from  (A. 32)  that 


T(p,q)L*  =  P*(p)q::Aj:A  (q+1)(s+1)  n  (A  -A. x (q+1)("i+1) 

P  P  i<p  P  1) 


*  .  11  (y-j-y,*)  n 

j >1>mp_i+q+1  J 


s+l  n.+l 

Mi  n  (y.-A.)  1 

J'>mp-l+q+1  J’>mp_l+q+1  J  1 

i<p-l 


• .  n  (u-x  )q+1 

J>mp_i+q+l  J  P 


(A. 33) 


Now  let  us  apply  the  operator  Tj|+1  2  (y  to  both  sides  of  (A. 33). 


It  follows  from  (A. 30)  that 


Vyq«  (VT(p>I,)  ■  T<p-o+1>- 


(A. 34) 


form-116  the  ri'9ht  Side>  Say’  R(p’q)  °f  (A'33)  in  the  snowing 

(A. 35) 


R(p,q)  ‘ 


where 


Rw  =  p*(p)q::Anq:An{q+1)(s+1)  n  (a  -a  )(q+1)K+1) 

i  p  p  v  p  Ai '  1 

i<P 


*  11  .  ,+q+2)  .  H  (y.-y,)yf+1 


j>m 


+n+?  J  Mm_  ,  +q+2;  .  .  n  +n+? 

_2+q+2  p-1  j>i>m  j+q+2  'l  ^  mp+l+q+^ 


"  ^  i<p,  (V1+q+rt)ni+1.  H  (Mj-X,)""’ 

1  <P~1  P-1  M  j>m^  -j+q+2  J  1 

i<p-l 


J  >rnp_  1  +q+2  J  i  <p-l  mp- 


.  *  (uj-xp)q+1 

J>rnp+i+q+2  J  p 


(A. 36) 
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Application  of  Tq+^  to  (A. 35)  yields  a  sum  of  terms  containing 

P-l  q 


derivatives  of  (ym  2-X  )q+1  with  respect  to  pm  2 

P-l  ^  ^  P-l 


of  orders 


0  up  to  q+1.  Thus,  we  can  write 


rq+l 


i+q+2R(p.q)  =  (qH):p^i+q+2  Rj 


+  ^mp_i+q+2"Xp^R2’ 


(A. 37) 


where  the  second  term  on  the  right  hand  side  is  the  sum  of  the  terms 

with  derivatives  of  (y  .  l0-A  )q+*  of  order  less  than  q+1.  There- 

mp-l  q  2  p 

fore  by  substituting  y  .  .0=A„  in  (A. 36)  and  in  (A. 37)  we  get 

mp_l+cl+^  P 


Tq+1 
mp_i+q+2  '“p 


(xn)R(p,q)  =  (q+1)1.  Apq+1  p*(p)q::xpq:xp(q+1)(s+1) 


n  (a  -x.)(q+1^ni+1)  n  (y.:-xp)  n  ^i"Pi)xp+1 

i<p-l  p  J>mp-l+q+2  J  P  j  >i  >nrip_  ^+q+2  J  p 


s+1 


n.+l 

n  n  ^  vxi ) 

j>rn  !+q+2  J  i<p-l  p 


n.+l 

n  (y,-A.) 
J*>mp-l+q+2  J 


n  (v<-xjq+1 


J>mp_i+Ci+2 


j  p 


or 


TVi+R*2  (Ap)R(p,q)=(c|+l)::x 


(q+1):  n  ( A„-A. ) 

i<p-l 


(q+2)(nf+l) 


•  n  (pf»  )qt2  p*(p)x  (q+2)(stl)  n  W.-X,)  1 

J>mp_l+Cl+2  J  P  "  j>mp_i+q+2  J 

i<p-l 


n.+l 
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n  (y,-x  )q+1  n  (y.-y-)  n  y. 
J>tnp-l+q+2  V  j>i>mp_1+q+2  J  j>mp_1+q+2  J 


By  comparing  this  with  (A. 33)  we  see  that 

T^+q+2  Up)R(p.q)  =  R(p,q+1).  (A. 38) 

In  view  of  (A. 34)  we  have 

T(p,q+1)  L*  =  R(p,q+1) .  (A. 39) 

* 

By  applying  the  operator  T  (j,p,q+l)  to  both  sides  of  (A. 39)  we  obtain 
(A. 26)  for  jng, • . . »np_j ,  and  np  replaced  by  q+1.  This  completes  the 

proof  of  the  lemma. 
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APPENDIX  B 


Details  of  the  procedure  for  numerical  differentiation  according 
to  the  method  of  this  report  will  be  described  in  a  forthcoming  report 
that  will  contain  examples  of  synthetic  data  as  well  as  applications  to 
experimental  records.  Here  we  include  only  a  few  illustrations  that 
compare  the  accuracy  of  this  method  with  the  accuracy  of  moving 
polynomial  arcs,  Butterworth  filter,  and  spline.  These  methods  were 
selected  for  comparison  because  they  are  frequently  applied  to  problems 
of  interest  to  BRL.  Furthermore,  the  corresponding  computer  programs 
are  available.  Also,  spline  approximation  yields  one  of  the  most 
accurate,  if  not  the  most  accurate,  derivatives  of  all  the  commonly 
used  methods  for  numerical  differentiation.  Comparison  with  Kalman 
filter  is  not  included  here  since  the  results  of  this  method  are  highly 
dependent  on  the  assumed  dynamics.  If  linear  equations  are  assumed 
that  are  satisfied  by  a  sine  function  of  a  given  frequency,  Kalman 
filter  produces  good  results  for  the  corresponding  synthetic  data. 
However,  if  the  same  dynamic  model  is  applied  to  synthetic  data  obtained 
from,  say,  an  exponential  function,  estimates  of  derivatives  are  poor. 
Illustration  of  these  cases  will  be  provided  in  the  forthcoming  report. 

The  synthetic  data  used  here  is  described  in  Table  1.  Here  column 
1  numbers  the  cases  from  1  to  6  for  convenience  of  reference.  The 
corresponding  functions  x(t)  are  specified  in  column  2.  The  last  case 
here  is  the  Bessel  function  of  the  first  kind  of  order  zero.  The 
values  of  the  functions  were  computed  at  points  in  the  intervals  listed 
in  column  4,  and  column  3  shows  the  step  size  for  selecting  the  points 
in  the  respective  intervals.  Pseudorandom  white  gaussian  noise  with 
zero  mean  and  the  standard  deviation  a  given  in  column  5  was  added  to 
each  value,  and  then  various  methods  for  numerical  differentiation  were 
applied  to  the  noisy  data. 

Table  1.  Synthetic  Data 


1 

2 

3 

4 

5 

Case 

x(t) 

A 

I 

a 

1 

sin  2irt 

.004 

[0,1] 

.02 

2 

sin  2irt 

.004 

[0,1] 

.05 

3 

e1 

.01 

[0,5] 

.01 

4 

4t(ltt+t2+t3) 

.004 

[0,2] 

.32 

5 

sin  2irt+.lsin  lOirt 

.004 

[0,1] 

.05 

6 

V‘> 

.01 

[1.6] 

.01 

51 


The  methods  used  are  listed  at  the  top  of  Table  2.  Here  moving 
polynomial  arc  corresponds  to  a  cubic  polynomial  fitted  to  either  11  or 
17  data  points  as  indicated.  The  derivatives  were  evaluated  at  the 
midpoint  of  this  span.  Thus,  derivatives  at  a  few  points  at  the 
beginning  and  the  end  of  the  data  sequence  are  not  available. 

The  Butterworth  filter  applied  here  corresponds  to  the  transfer 

function  tt/(s  +2tts  +2tt  s+ir  ).  This  program  does  not  provide  first 
derivatives  at  65  data  points  at  the  end  of  the  data  sequence  when 
A  =  .01  and  at  163  data  points  when  A  =  .004.  Additional  points  are 
lost  when  higher  derivatives  are  calculated.  The  method  of  this  report 
provides  derivatives  at  every  data  point  with  appropriate  values  of  s 
in  (3.2)  and  in  the  corresponding  expressions  for  higher  derivatives. 

The  bulk  of  the  derivatives  are  computed  at  the  midpoint  of  the  span 

of  formula  (3.2),  i.e.,  for  s  =  -  jj^-J  .  Estimates  of  the  bounds  as 
described  in  Section  5  also  depend  on  s.  In  the  examples  of  this  appen¬ 
dix  we  use  only  s  =  i.e.,  the  RMSE's  of  the  numerical  derivatives 

are  obtained  here  by  comparing  the  derivatives  computable  by  (3.2)  with 
this  value  of  s  against  the  exact  derivatives.  Also  an  estimate  of  the 
bound  for  the  standard  deviation  of  the  error  is  obtained  by  (5.13)  for 
only  these  points. 

The  double  rows  labeled  1-6  in  Table  2  show  the  percentage  errors 
corresponding  to  the  data  sets  1-6  of  Table  1  obtained  by  the  methods 
shown  at  the  top  of  Table  2.  The  RMS  of  the  differences  between  exact 
and  approximate  first  derivaties  is  divided  by  the  RMS  of  the  exact 
derivative  and  the  ratio  is  multiplied  by  100.  This  yields  the 
percentage  error  of  the  first  derivative  for  the  respective  method. 

This  error  is  shown  in  the  first  row  of  a  pair  labeled  x'  .  The  second 
row  of  each  pair  labeled  x/;  contains  the  percentage  errors  of  the  second 
derivative  obtained  in  the  same  way.  This  table  shows  that  the  method 
of  the  report  (current  method)  when  applied  to  these  synthetic  data 
yields  far  more  accurate  derivatives  than  any  of  the  other  methods. 

For  instance,  in  cases  1-3  the  first  derivative  obtained  by  our  method 
is  from  two  to  more  than  four  times  as  accurate  as  by  the  second  best 
method  (spline),  while  the  error  in  the  second  derivatives  is  lower  by 
more  than  an  order  of  magnitude.  In  the  remaining  cases  the  current 
method  also  produces  much  better  results  than  other  methods  except  for 
case  6  where  derivatives  obtained  by  the  Butterworth  filter  are  more 
accurate . 


Bl 
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Table  2.  Percentage  Error 


Method 

Cas^v 

Moving  Polyn. 

Arc,  11  pts 

Moving  Polyn. 

Arc,  17  pts. 

Butterworth 

Fi 1  ter 

Spl ine 

Current 

Method 

X' 

1 

X" 

25 

14 

88 

6 

1.3 

275 

99 

98 

48 

1.4 

no 

X  X 

^  N 

62 

33 

88 

15 

3.4 

688 

248 

98 

120 

3.5 

x' 

3 

X" 

.46 

.25 

3.4 

.010 

.004 

13 

4.7 

6.9 

2.2 

.004 

X' 

4 

X" 

9 

5 

16 

2.1 

.23 

259 

93 

30 

45 

.80 

X/ 

5 

55 

30 

91 

13 

8 

X" 

89 

255 

99 

46 

18 

X' 

6 

57 

31 

5.8 

13 

7.8 

X" 

2341 

- i. 

8384 

8.4 

403 

33 

Table  3  illustrates  how  an  approximation  of  the  bounds  of  the 
error  compares  with  the  actual  RMSE's.  Here  the  first  row  lists  actual 
percentage  errors  of  the  first  and  second  derivatives  for  the  various 
cases  listed  in  the  headings  of  the  column-pairs.  This  row  is  the  same 
as  the  last  column  in  Table  2.  The  second  row  here  contains  the  ratios 
of  the  bounds  cr(r)  (r=l,2)  computed  by  (5.13)  to  the  RMS  of  the  respec¬ 
tive  exact  values  converted  to  percent.  We  see  that  (5.13)  yields 
bounds  about  2  to  10  times  higher  than  the  respective  RMSE's  except  for 
the  second  derivative  of  case  6  (Bessel  function)  where  the  estimate  of 
the  bound  is  much  too  low. 


Table  3.  Error  Bounds 


1 

2 

3 

4 

5 

6 

X7 

X77 

X7 

X77 

X7 

X77 

X7 

X77 

X7 

X77 

X7 

X77 

1.3 

1.4 

3.4 

3.5 

.004 

.004 

.23 

.80 

8 

18 

7.8 

33 

3.5 

2.4 

8.9 

6.2 

.02 

.02 

2.1 

3.0 

29 

91 

10 

4.5 
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